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ABSTRACT 

We analyse the concentration of solid particles in vortices created and sustained by radial buoy¬ 
ancy in protoplanetary disks, i.e. baroclinic vortex growth. Besides the gas drag acting on particles 
we also allow for back-reaction from dust onto the gas. This becomes important when the local dust- 
to-gas ratio approaches unity. In our 2D, local, shearing sheet simulations we see high concentrations 
of grains inside the vortices for a broad range of Stokes numbers, St. An initial dust-to-gas ratio of 
1:100 can easily be reversed to 100:1 for St = 1. The increased dust-to-gas ratio triggers the streaming 
instability, thus counter-intuitively limiting the maximal achievable overdensities. We find that par¬ 
ticle trapping inside vortices opens the possibility for gravity-assisted planetesimal formation even 
for small particles (St = 0.01) and low initial dust-to-gas ratios (1:10 4 ). 

Subject headings: planets and satellites: formation, protoplanetary disks, accretion disks - hydrody¬ 
namics - instabilities 


1. INTRODUCTION 

A quantitative prescription of planetesimal formation 
is one of the key issues of planet formation theory. Mod¬ 
els of simple collisional sticking are controversially dis¬ 
cussed, both conceptually and also in the explanation 
of current properties of asteroid and Kuiper belt obejcts 
(Weidenschilling 2011; Morbidelli et al. 2009; Nesvorny 
et al. 2011; Bottke et al. 2005). Whether small dust grains 
stick to one another, bounce or fragment depends on 
their size and their relative velocities. In general, theory 
predicts that collisional velocities rise as particles grow, 
which holds for particles with Stokes number smaller 
than 1. The Stokes number is defined as St = Dtj, where 
Q is the Keplerian frequency, and Tf is the stropping 
time of a particle; both drift velocity and turbulence- 
induced relative velocities have a maximum for St=l, 
i.e., when the stopping time is of the order of an or¬ 
bital period (see e.g. review by Weidenschilling & Cuzzi 
1993). At 5 AU, St=l corresponds to meter-size objects. 
Fragmentation occurs at velocities of only a few ms -1 , 
which limits particle sizes to mm-cm (Wurm & Blum 
2000; Brauer et al. 2008; Birnstiel et al. 2010). Guttler 
et al. (2010) and Zsom et al. (2010) introduced another 
boundary, the so called bouncing barrier where parti¬ 
cles hit one another and bounce without mass transfer. 
At even smaller size scales, Okuzumi et al. (2011a,b) 
found the charge barrier, where small particles are pre¬ 
vented from approaching one another due to the elec¬ 
tric charges built up through collisions with free elec¬ 
trons. Birnstiel et al. (2012) determined the sizes that 
particles can obtain locally as their growth is limited by 
radial drift and collisional destructions. They find al¬ 
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most independently on distance from the star a max¬ 
imum Stokes Number of about 0.01 - 0.1 which corre¬ 
sponds to 1 -10 cm at 1AU, 0.3 - 3 cm at 10 AU and 
0.3 - 3 mm at 100 AU (see figs. 5 & 6 in Birnstiel et al. 
2012 ). 

To form large planetesimals these difficulties need to 
be circumvented. One proposed method is gravitational 
instability (Safronov 1972; Goldreich & Ward 1973; Jo¬ 
hansen et al. 2006a, 2007): When a sufficient amount of 
particles is close enough together, their mutual attrac¬ 
tion can trigger gravitational collapse, rapidly forming 
large planetesimals that then sweep up small particles 
(Lambrechts & Johansen 2012). Different methods to 
capture dust have been studied, such as zonal flows by 
Johansen et al. (2011) and Dittrich et al. (2013), pressure 
rings around stars (Whipple 1972; Klahr & Lin 2001, 
2005), convection cells (Klahr & Henning 1997) or vor¬ 
tices either numerically (Tanga et al. 1996; Johansen et al. 
2004; Lyra et al. 2008a, 2009a,b; Meheut et al. 2012a,b) 
or analytically (Barge & Sommeria 1995; Chavanis 2000; 
Klahr & Bodenheimer 2006; Chang & Oishi 2010; Lyra & 
Lin 2013). 

From these studies it has become clear that dust can 
easily concentrate in anti-cyclonic vortices. However, 
feedback of dust on the vortical flow, which can poten¬ 
tially frustrate particle concentration, has not been stud¬ 
ied in detail. This is the topic of the present study. 

Besides Lyra et al. (2008a, 2009a,b); Meheut et al. 
(2012b); Ataiee et al. (2013); Zhu & Stone (2014), and 
Zhu et al. (2014), who analyzed particle trapping in vor¬ 
tices excited by the Rossby wave instability, other works 
did not choose a particular vortex formation mecha¬ 
nism. In our study, vortices are naturally produced by 
the radial stratification of disks (Klahr & Bodenheimer 
2003; Petersen et al. 2007a,b; Lesur & Papaloizou 2010; 
Lyra & Klahr 2011; Raettig et al. 2013). The growth of 
the vortices occurs proportionally to the radial buoy¬ 
ancy frequency (aka Brunt-Vaisala frequency) squared 
N 2 (Raettig et al. 2013) and the local thermal relaxation 
time scale as argued by Lesur & Papaloizou (2010), nu¬ 
merically confirmed by Raettig et al. (2013) . The buoy- 
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ancy frequency, which is a function of the relative scale 
height H/r, logarithmic radial entropy and pressure 
gradients jS k,Pp and the adiabatic index 7 : 

A linear theory to explain this behaviour has been put 
forward by Klahr & Hubbard (2014) who identify the 
instability as a radial "convective overstability" in accre¬ 
tion disks. Lyra (2014) performed 3D simulations of this 
"convective overstability" and showed that in the non¬ 
linear phase indeed vortices were emerging from the 
flow. 

In this paper we set out to assess how efficiently 
dust can be concentrated in vortices enforced by real¬ 
istic values for the radial stratification in temperature 
and density, also with a plausible value for the ther¬ 
mal relaxation time. We do this as a first step via two- 
dimensional, local simulations. Ultimately, only three- 
dimensional stratified runs will be able to include all 
relevant physics from vortex stability (Barranco & Mar¬ 
cus 2005; Lesur & Papaloizou 2009; Lyra & Klahr 2011) 
to sedimentation of dust. Nevertheless, here we find 
even in 2D that particle concentration in the vortices 
is sufficiently intense that back-reaction of the particles 
onto the gas motion has to be considered. We find that 
the streaming instability (Youdin & Goodman 2005; Jo¬ 
hansen et al. 2006a) is triggered, counter-intuitively lim¬ 
iting the maximum dust-to-gas ratio, and severely per¬ 
turbing the gas flow inside the vortex. In some cases 
we see vortices getting disrupted and later reforming, 
starting a new cycle of particle concentration. 

In simulations similar to ours but three-dimensional 
and more numerically expensive, where particles were 
trapped in zonal flows of magneto-hydrodynamical ori¬ 
gin and including self-gravity, the overdensities that we 
found in the present paper did already lead to gravita¬ 
tional collapse (Johansen et al. 2006b, 2007) and to the 
formation of planetesimals. On the other hand 2D sim¬ 
ulations including particle feedback are ideally suited to 
study a wide parameter range to learn whether stream¬ 
ing instability can be triggered and how it affects vortex 
stability. 

The paper is structured as follows. We first review the 
underlying physics of dust-motion in a gas disk, specif¬ 
ically in a vortex including dust-gas interactions. Then 
we describe the numerical setup. The results of our sim¬ 
ulations are given in Section 4. Here we especially look 
at the reached dust-to-gas ratios compared to the aver¬ 
age initial dust-to-gas ratio. In Section 5 we discuss the 
prospects to form planetesimals via gravitational frag¬ 
mentation of the dust enhancements in vortices and in 
Section 6 we present our measurements of collisional 
speeds among dust grains. Finally, we summarize and 
conclude in Section 7. 

2 . PHYSICAL BACKGROUND 

The evolution of the gas component of the disk in the 
Pencil Code is given by the Navier Stokes Equation con¬ 
taining stellar gravity as well as the virtual forces of 
the rotating and shearing box as explicit terms (Lyra & 
Klahr 2011). Dust grains on the other hand are evolved 


by solving the equations of motions for Lagrangian par¬ 
ticles (see Section 3). 

In general the gas and dust feel the same external 
forces except the pressure force f v = -p " 1 Vp, where p g 
is the gas density and p the pressure. This term only 
affects the gas. For instance the global pressure gradi¬ 
ent in the gas leads to a sub-Keplerian orbital gas veloc¬ 
ity u .The corresponding buoyancy force for a particle is 
fp,s = “Ps 1 ^? 7 / where p s is the material density of the 
solid material, can be neglected because p s p g . Since 
the particles do not feel this global pressure, they need 
to orbit at Keplerian velocities in order to be in centrifu¬ 
gal balance with stellar gravity. The resulting veloc¬ 
ity difference between gas and large particles acts as a 
headwind on the particles which decreases their angu¬ 
lar momentum leading to radial drift inward. Smaller 
particles get dragged along by the gas and therefor feel a 
net acceleration towards the star also making them drift 
inward. 

The time on which the dust particles adjust to the gas 
velocity is the friction time. For subsonic Epstein drag 
this is r s = pss(pgC s ) -1 (Weidenschilling 1977) which 
depends on particle size s and local sound speed c s . 
The Epstein regime considers particles smaller than the 
mean free path of gas. Larger particles have to be 
treated in the Stokes regime. But as long as relative ve¬ 
locities between dust and gas are small, it is possible 
to describe the coupling force to be linear with respect 
to the relative velocity f = m^r independent on the de¬ 
tailed calculation of r s . Especially it is not necessary to 
determine the shape and density of the particles or the 
gas density and sound speed until one asks for the par¬ 
ticle size corresponding to a given friction time. 

This particle gas coupling is usually expressed in 
terms of the dimensionless Stokes number St = Qt s . 
Particles of different size, shape and density, but ulti¬ 
mately with the same St will behave the same as far as 
aerodynamics are concerned. 

The friction time and thus the Stokes number is a 
function of the local gas density. We nevertheless as¬ 
sume the Stokes number to be constant, which is justi¬ 
fied as the gas density fluctuates only by a few percent 
in the sub-sonic turbulence we consider in our simula¬ 
tions. 

Since our simulations are only two-dimensional we 
can not consider vertical settling of particles. The fact 
that the real 3D gas disk is vertically stratified and thus 
particles of a given size have a vertically varying friction 
time, has to be discussed in more detail. One can inter¬ 
pret this approximation as focussing on a narrow verti¬ 
cal range around the midplane, in which neither the gas 
nor the particle stratification are too strong. Then one 
can assume that also the friction time and the dust to 
gas ratio has no strong vertical variations. This assump¬ 
tion is a little less stringent than saying that all particles 
are actually in the midplane. 

The general behavior of particle drift in a non-laminar 
yet geostrophic flow 5 can be understood easily from a 
a simplified treatment of the Lagrangian equations of 

5 In a geostrophic flow the pressure gradients are in equilibrium 
with the centrifugal and Coriolis acceleration. 
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motion for gas and dust: 


TABLE 1 

Simulation setup and e max 


d t u = F-—Vp--^-(u-v) (2) 

Ps Ps T * 

dtv = F- — (v-u), (3) 


where v is the particle velocity and F collects the terms 
that are the same for both dust and gas, e.g., gravita¬ 
tional force from the central star. The second term of 
the equations describe the pressure force, and the third 
term is the friction force between gas and dust particles. 
With some math it can be shown that by subtracting Eq. 
(2) and (3) and assuming \dtu-dtv\ <C \u-v\/r s par¬ 
ticles will move towards regions with higher pressure 
v = u + t s V p even if the pressure maximum is tiny and 
the profile relatively smooth (see e.g. Klahr & Lin 2001). 

The vortices we consider are anti-cyclonic vortices 
and have a slightly lower epicyclic frequency than the 
Keplerian orbital frequency. Since there is a geostrophic 
balance between the Coriolis and the pressure force a 
high pressure region inside the vortex is created. Parti¬ 
cle accumulation inside a vortex then basically works 
the same way as the radial drift in an accretion disk 
works. The particles do not feel the pressure gradient 
inside the vortex and therefore their epicyclic frequency 
equals the Keplerian frequency. Yet, the rotation fre¬ 
quency of a pressure supported vortex is smaller than 
the epicyclic frequency. Thus the headwind from the 
gas causes the particles to lose eccentricity and forces 
the particles to spiral towards the center of the vortex. 
For an in depth analysis of vortex capturing including a 
comparison of all forces acting on the gas and dust in¬ 
side the vortex see Adams & Watkins (1995); Barge & 
Sommeria (1995); Tanga et al. (1996); Chavanis (2000) 
and Johansen et al. (2004). 

If the dust density becomes comparable to the gas 
density the drag forces from particles onto the gas can 
no longer be neglected as long as the St is smaller than 
« 10. This back-reaction can alter the motion of the 
gas and also leads to even higher dust concentrations 
through the streaming instability (Youdin & Goodman 
2005). The last term of equation (2) represents the back- 
reaction from dust grains onto the gas. 


3. NUMERICAL SETUP 

We perform two dimensional shearing sheet simula¬ 
tions with the Pencil Code where we consider the 
vertically integrated densities Z instead of the three- 
dimensional densities p. The Euler-equations for the gas 
are solved on a Cartesian grid, identical to the setup 
of (Raettig et al. 2013), but now augmented by a term 
for the particle feedback on the gas e (u - v) /r s with the 


run 

St 

^d,o/Z g/ o 

feedback 

£max 

Fraction of 

with £ > 0.01(%) 

Fraction of 

with £ > 1(%) 

NF1 

0.01 

1:100 

no 

6.26 

86.64 

2.88 

NF2 

0.05 

1:100 

no 

636.00 

92.31 

75.63 

NF3 

1 

1:100 

no 

843.70 

79.71 

1.26 

NF4 

20 

1:100 

no 

46.81 

99.90 

99.80 

FI 

0.01 

1:100 

yes 

1.07 

70.40 

0 

F2 

0.05 

1:100 

yes 

3.86 

95.24 

2.13 

F3 

1 

1:100 

yes 

77.33 

98.73 

83.74 

F4 

20 

1:100 

yes 

0.73 

79.82 

8.09 

DG1 

0.05 

1:1000 

yes 

1.15 

86.74 

0 

DG2 

0.05 

1:10000 

yes 

0.70 

78.30 

0 

DG3 

1 

1:1000 

yes 

11.53 

98.46 

56.37 

DG4 

1 

1:10000 

yes 

4.17 

97.30 

76.15 


TABLE 2 

Collision al velocities 

St 

£ o 

DC 

uw~ 2 ) 

^coll 

Cs 

(/1(T 3 ) 

Urms 

V0LC S 

Prms 

y/oLCs 

^coll 

yfocc s 

(/1(T 2 ) 

A ^coll 

y/lXC s 

(/1(T 2 ) 

u coll,OC 

y /* c s 

1 

10 “ 2 

1.03 

2.77 

2.84 

1.37 

2.73 

2.56 

1 

1 

10" 3 

1.56 

3.17 

4.36 

3.61 

2.54 

3.34 

1 

1 

10 -4 

1.54 

1.25 

3.88 

2.27 

1.01 

4.01 

1 

0.05 

10" 2 

1.56 

2.86 

3.27 

2.41 

2.29 

0.59 

0.31 

0.05 

10“ 3 

1.57 

2.71 

4.29 

3.80 

2.16 

1.28 

0.31 

0.05 

10- 4 

1.54 

2.08 

3.88 

3.04 

1.68 

0.09 

0.31 


dust-to-gas ratio e = Z^/Zg 
DZ 

' + (u • V)Z g = -ZgV • m +/d( Z g ), 

^ + (m- V)ii = -^-Vp- 2 U 0 (z x 1 1 ) 

+ y^0 u xV + ^ x + /v(w,Z) 


( 4 ) 


Z g T s 


Ro 
(u-v), 


( 5 ) 


Vs 

—- + (u • V)s = _ 
Vt v J Z 


(T-To) 


+ 


yV.(KVT)-^ Tcooi 

^ } + A( S ). ( 6 ) 


r 0 (7-I) 


Here, Z g is the gas density, u is the deviation of the 
gas velocity from the Keplerian value, s the entropy, T 
the temperature, c v the specific heat at constant volume, 
Tcool is the thermal relaxation time scale, and K the heat 
conductivity. The symbol 


V 

Vt 



(0)d_ 
y dy 


( 7 ) 


represents the Keplerian derivative where u ^ = 
-3/2i7 0 x and Uq is the Keplerian frequency at the or¬ 
bit of the shearing box Rq. The radial deviation from 
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Rq is given by x and the linearized azimuthal direc¬ 
tion is now measured in terms of y. Boundary condi¬ 
tions are periodic in y and shear-periodic in x. Fur¬ 
ther terms in the equations are diffusion terms to en¬ 
sure numerical stability of the finite difference Pencil 
Code /D(Sg) / /v(w,S) / /x(s) and radial stratification of 
the disk in terms of pressure and entropy for radially 

constant density: /5 = The stratifica¬ 

tion term (5 occurs in the linearized component of radial 
buoyancy (Eq. (5)) and in the term for radial transport 
of entropy (Eq. (6)) as derived in Lyra & Klahr (2011). 
These terms in combination with thermal relaxation are 
the driver of vortices as found in Petersen et al. (2007a). 
As described in their paper gas that is moving radially 
outward is generally warmer and of lower density than 
gas moving inward. The density mismatch across the 
vortex leads to a mismatch in buoyancy and thus a vor¬ 
tex will feel a torque accelerating it. The maximal ampli¬ 
fication occurs when thermal relaxation is on the same 
order as the internal rotation period of the vortex. 

The dust grains are modeled with a particle approach. 
For each individual particle we solve the equation of 
motion including gas drag. We do not allow for self¬ 
gravity so far. In principle, for an individual particle it 
does not matter if there are other particles in the simula¬ 
tion or not. However, in the simulations where we allow 
for back-reaction of the particles on the gas, there is an 
indirect influence of one particle on the other particles 
due to the altered gas velocity. 

The gas disk is co-rotating with Keplerian velocity at 
the co-rotational radius Tq. But for simulations involv¬ 
ing both dust and gas we need to include a velocity off¬ 
set due to the pressure gradient. The pressure gradient 
is balanced by the Coriolis force which is the linearized 
expression for the radial centrifugal acceleration in the 
local co-rotating system 


2 f2u y = 


1 dp Hdlnp 

V O C C TCj . 

£ g or r dlnr 


( 8 ) 


The deviation by which the gas velocity is lower than 
the Keplerian velocity u k = f2r is (Nakagawa et al. 1986) 




2 din p 
din r' 


( 9 ) 


which leads to a sub-Keplerian gas velocity and a devi¬ 
ation from the Keplerian velocity u k of 


Uy - 0.5 



2 


din p 
din r 


Qr = y w K . 


( 10 ) 


This velocity deviation y is added to the simulations ar¬ 
tificially. 

For each particle i the equation of motion is solved 
individually via 


dv$ 

dt 


= 2f2v*yX— -u (*^)) , (H) 


position of the dust particle. The index * is omitted from 
now on. 

For our 2D simulations we choose y = -0.01 which 
corresponds to a pressure gradient of /3 = 2 and a disk 
aspect ratio of h = 0.1. The terminal velocities of the 
particles for a given Stokes number in a laminar disk 
are now in terms of y (Weidenschilling 1977) 


v 


x — 


Vy = 


k 

st + sr 1 ' 

st 2 ’ 


( 12 ) 

(13) 


These are the best assumptions one can do for the initial 
velocities of particles in our simulations. 

Our physical domain spans 4 disk scale-heights, ±2H, 
around the co-rotational radius in radial-direction (x- 
axis) and 16H in azimuthal-direction (y-axis). The grid 
itself consists of 288 2 grid cells 6 . We choose jS = 2 for 
the entropy and pressure gradients which is a relatively 
strong gradient compared to gradients we expect in pro- 
toplanetary accretion disks are between = 0.5 and 
j6 = 1. However, Raettig et al. (2013) found that the 
general behavior of vortices is the same for weak and 
strong entropy gradients. The development of the vor¬ 
tices is merely faster with stronger gradients and we can 
scan a large parameter space with a reasonable amount 
of computing time. For a first estimation the strong gra¬ 
dient is sufficient, but in our future 3D studies including 
self-gravity we will use (5 = 1 because then we aim for 
quantitative results on planetesimal formation rates and 
mass spectrum. 

We first evolve the disk for 200 local orbits without 
particles. This way we make sure that a fully grown, 
long lived vortex has developed before we put in par¬ 
ticles. Physically, this corresponds to various scenarios: 

A) a vortex developing in the outer parts of the disk and 
then migrating inwards through regions with particles, 

B) the growth of particles into a size regime where they 
get trapped by the nearest vortex or C) a radial flux of 
particles close to bouncing and / or drift barrier that en¬ 
ters the vicinity of the vortex. 

After the initial 200 orbits, we distribute 400 000 parti¬ 
cles randomly in the disk. This corresponds to 4-5 par¬ 
ticles per grid cell initially. That way we minimize nu¬ 
merical effects that can arise if there are not enough par¬ 
ticles in the computational domain, e.g. the effect a sin¬ 
gle particle has can be extremely overestimated if not 
enough particles are considered. Whenever we refer to 
times in this paper, we mean time elapsed since particles 
were put into the simulation. At the time when the dust 
to gas ratio approaches unity and streaming instability 
is triggered, 500 or more particles are in one grid cell. 
Following Johansen & Youdin (2007) these numbers are 
safely beyond the minimum particle number per grid 
cell needed to achieve numerical convergence. 

Starting our simulations with particles included 
would have little changed our results as the initial and 
average dust to gas ratios are much too low in order for 
the dust to have impact on the gas dynamics right from 


where = v - y is the particle velocity corrected by 
the velocity offset and u it the gas velocity at the 


6 The number of 288 grid cells comes from the architecture of our 
computational cluster which needs a multiple of 12 to parallelize the 
Pencil Code efficiently. 
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Fig. 1.— Maximum dust-to-gas ratio, £ max , for simulations without 
particle feedback. The different lines represent the different particle 
sizes: dotted (black) line: St = 20, dash-dotted (green) line: St = 1, 
solid (blue) line: St = 0.05, and dashed (red) line: St = 0.01. Almost 
all particles of intermediate size (St = 1 and St = 0.05) accumulate in 
the vortex. Strongly coupled particles (St = 0.01) accumulate only par¬ 
tially, because they also couple to the gas outside of the vortex. Large 
particles (St = 20) hardly couple to the gas at all and therefore are not 
affected by the vortical motion. 

the beginning. Retrospectively this assumption was 
later justified in simulations (see Section 4.2) in which 
the vortex gets disrupted via the streaming instability 
and then forms again, this time in the presence of dust 
grains. This second growth phase does not differ from 
the initial growth phase without dust grains present. 

Each of the 400000 particles represents one super¬ 
particle, a collection of several particles, all of the same 
size, and a given mass according to the initial dust-to¬ 
gas ratio. 

Generally we have an initial dust-to-gas ratio £q of 
1:100 which means that the disk consists out of 1% solid 
material and 99% gas. Note that in two dimensions the 
Pencil Code assumes surface densities instead of vol¬ 
ume densities. Therefore the dust-to-gas ratios we talk 
about in this paper refer to e = Z^/Eg rather than to 
e = Pd/Pg- To simulate different particle sizes we use 
different Stokes numbers of St = 0.01, St = 0.05, St = 1 
and St = 20. By only defining the Stokes Number we 
are able to cover particles both in the Epstein as well 
as in the Stokes regime as friction forces are linear with 
respect to velocity. At 5 AU this corresponds to par¬ 
ticles between 3mm and 6m for Z g = 300gem -2 and 

p s = 2gem -3 . The parameters for all our simulations 
can be seen in Table 1. A separate simulation is carried 
out for each St. 

Particles growth due to collisions is not treated in our 
simulation. Physically particles would have a growth 
time of several thousand years, therefore we neglect this 
effect for the sake of keeping the simulation simple and 
easy to be evaluated. 

4. RESULTS 
4.1. No Particle Feedback 

We first consider simulations where we do not include 
back-reaction from the particles, effectively setting p ^ = 
0 in Equation (5). That means that gas drags particles 
along, but even high dust densities do not affect the gas 
velocity in any way. The purpose of these simulations 
is to serve as comparison to the later models in which 



Fig. 2.— Azimuthally averaged gas velocity and particle density for 
St=20 particles. The particles accumulate in the zonal flow structure. 

feedback is included. 

Figure 1 shows the maximum dust-to-gas ratio for 
these simulations. We see from the almost constant 
dash-dotted green line that all of the available solid ma¬ 
terial in the St = 1 particle simulations accumulates in 
the vortex. This corresponds to a particle concentration 
for St=l particles of c « Z^/Z^o ~ 12000 and £ ~ 320. 
St = 0.05 eventually also all accumulate within the vor¬ 
tex and then are kept trapped there. The maximum par¬ 
ticle concentration is reached sooner for St = 1 particles 
than for St = 0.05 particles, which can be attributed to 
their higher radial drift velocity. Even for St = 0.01 par¬ 
ticles the dust-to-gas ratio increases to £ « 1 . 

We see two different behaviors. For St = 0.05 and 
St = 1, basically all particles are collected within the vor¬ 
tex, although St = 0.05 take 10 times longer, yet in both 
cases concentrations of about 4 orders of magnitude are 
reached. For St = 0.01 concentration inside the vortex is 
only by 2 orders of magnitude and St = 20 particles con¬ 
centrate by one order of magnitude, but interestingly 
outside the vortex in correlation with the zonal flow 
related to the radial position of the vortex (see Fig. 2). 
These zonal flows show up in the azimuthally averaged 
rotation profile, both in magneto hydrodynamical sim¬ 
ulations (Lyra et al. 2008b; Johansen et al. 2009; Dittrich 
et al. 2013) and also in our hydrodynamical simulations 
of vortices. The zonal flows are deviations from the sub- 
Keplerian mean rotation profile confined by long-lived 
radial pressure variations (geostrophic balance). 

The different behavior of particle accumulation can 
be explained through the gas coupling. Particles with 
St = 0.01 are strongly coupled to the gas. As an effect 
it takes much longer to capture the particles into the 
vortex and as a second effect the deviation between the 
actual vortex and a analytical solution like for instance 
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Fig. 3.— Each of the four panels show vorticity co z (in units of the local Keplerian frequency Qq), and dust-to-gas ratio e. The left-hand-side 
panels do not include particle feedback, whereas the ones in the right-hand-side do. The particle sizes correpond to Stokes numbers St — 20 (top) 
and St = 1 (bottom). 
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Fig. 4.— Same as Fig. 3, but for particles of St = 0.05 (top) and St = 0.01 (bottom). 
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Fig. 5.— Gas vorticity (left panel) and dust-to-gas ratio (right panel) 
for a simulation with St = 0.01 particles and back-reaction onto the 
gas. The elliptical vortical gas flow is distinguishable in the vorticity 
plot. There are strong accumulations of particles within the vortex. 
Although many particles spread out over the entire vortex, most par¬ 
ticles concentrate in the center of the vortex. The positive vorticity 
values in the vortex (light areas in the left plot) show the effect that 
particles have on the gas. Where the gas encounters obstacles, namely 
high particle concentrations, steep vorticity gradients develop. 


a Kida vortex (see discussion in Lyra & Lin 2013) pre¬ 
vents the particles from accumulating too densely. An 
additional effect to prevent too strong concentrations in 
3D vortices might be turbulence, triggered for instance 
by elliptical instability (Lyra & Lin 2013) and internal 
circulation (Meheut et al. 2012a). The St = 20 particles 
are the other extreme: they are hardly coupled to the 
gas at all, and therefore are only weakly affected by the 
vortical motion of the gas. In a stationary analytic vor¬ 
tex these particles would get trapped if one waits long 
enough. Yet, our numerical vortex is dynamically ac¬ 
tive in changing its strength and shape and additionally 
drifting in the radial direction. This vortex dynamics 
gives the large particles the chance to escape from the 
vortex again, but still stopped from radial drift in the 
zonal flow. 

4.2. Including Particle Feedback 


In the last Section we showed that there are significant 
particle overdensities within the vortex. This means that 
we cannot neglect the momentum transfer from dust 
onto gas. Because of this we now include dust back- 
reaction on the gas, using the full Equation (5). 

For the St = 20 particles there is no change with re¬ 
spect to the simulations without particle feedback be¬ 
cause overdensities stay below e = 1 (compare the left 
and right side of the top panel in Fig. 3). For these large 
particles the streaming instability can not be triggered 
at the present dust-to-gas ratio. All other tested parti¬ 
cles still accumulate in the vortex, yet differently: 

Particles with St = 1 concentrate more locally (in a 
smaller area) than smaller particles (bottom panel of 
Fig. 3 and Fig. 4). St=l particles are getting concentrated 
on time scales which are not longer than the dynami¬ 
cal time scale of the vortex (Barge & Sommeria 1995), 
thus they can easily follow the changes of the vortex 
in shape, strength and location. Smaller particles have 
much longer time scales to spiral into the center of the 
vortex - longer than the dynamical time scale of the vor¬ 
tex. As the attractor inside the vortex is changing its lo¬ 
cation on dynamical time scales (change of amplitude, 
shape and location of the vortex), the particles have no 
chance to ever catch up or follow it. 

An effect that occurs when e approaches unity is the 
streaming instability (Youdin & Goodman 2005; Youdin 
& Johansen 2007). Once particles concentrate, their lo¬ 
cally increased dust-to-gas ratio leads to a slower radial 
drift (Nakagawa et al. 1986). This produces a further en¬ 
hancement of solids since faster particles from slightly 
larger radii bump into the accumulation, like a traffic 
jam, which eventually results in streaming dust struc¬ 
tures. This is also the case in our simulations, although 
with our resolution it is possible to resolve only a part of 
the unstable wavelengths . Therefore we are not able to 
study the linear and nonlinear behavior of the streaming 
instability in all detail and model the correct growth- 
rates of the instability Yet, for the wavelength we re¬ 
solve we find the right threshold dust-to-gas ratio to 
trigger the instability See also the appendix for addi¬ 
tional resolution tests on the streaming instability 

Figure 5 shows the vertical gas vorticity co z in units 
of the local Keplerian frequency and dust-to-gas ra¬ 
tio e = Egfor a simulation with St = 0.01 particles 
after 300 local orbits (run FI). It is clear that the parti¬ 
cles accumulate inside of the vortex and follow the vor¬ 
tical motion. Where the concentration is highest they 
create strong maxima in the gas vorticity This is an ef¬ 
fect of the cyclonic rotation of the particle clumps which 
was already reported in earlier work (Lambrechts & 
Johansen 2012) and is a consequence of the conserva¬ 
tion of angular momentum under contraction which can 
turn anticyclonic motion in a rotating system to cyclonic 
motion. 

In general due to the back-reaction of the dust on the 
gas and the resulting streaming instability the initially 
elliptical gas streamlines are bent into more complex 
motions than in cases without back-reaction. 

In case of the St = 1 particles where the local parti¬ 
cle concentration, and therefore the back-reaction, are 
strongest, the vortex structure is disrupted. Because 
of this, the particle trapping mechanisms lose strength: 
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t= 99.37 


FIG. 6.— Vorticity co z in units of the local Keplerian frequency !?o (1st and 3rd column) and dust-to-gas ratio (2nd and 4th column) for St = 1 
particles for different points in time. At the first snapshot, the vortex is still clearly distinguishable. In the second snapshot, it has been disrupted 
strongly by the particle accumulation. As this particle accumulation spreads out, the vortex can slowly regain its shape (3rd snapshot) and form a 
large, yet still perturbed, vortex again (4th snapshot). 








































10 


Raettig et al. 



FIG. 7.— Maximum dust-to-gas ratio, £ max , for simulations with par¬ 
ticle feedback. The different lines represent the different particle sizes: 
dotted (black) line: St = 20, dash-dotted (green) line: St = 1, solid 
(blue) line: St = 0.05, and dashed (red) line: St = 0.01. Particles of 
St = 1 reach the highest dust enhancements and concentrate in a very 
local area in the vortex, while smaller particles spread out over the 
entire vortex. Therefore the overdensities reached are lower. 



Fig. 8.— Fraction of the entire dust mass in a specific dust-to-gas ra¬ 
tio or higher, e > 1 are reached for all particle sizes. For St = 1 particles 
more than 80% of all particles have concentrated in areas with larger e 
than 1. 



Fig. 9.— Fraction of the entire dust mass in a specific dust-to-gas 
ratio or higher, e > 1 are reached for St=0.05 particles with (dashed 
black line) and without feedback (solid blue line). 

pressure gradients across the vortex become shallower, 
Coriolis forces in the vortex grow weaker, and the parti¬ 
cles begin to escape the vortex. Because the local parti¬ 
cle concentration decreases, the large vorticity gradients 
flatten out again. This eventually leads to a new amplifi¬ 
cation of the vortex due to the background stratification 
and the process repeats itself (see series of snapshots in 



0 100 200 300 400 

time [27 c/£2 () ] 


Fig. 10.— Particle concentration (solid lines) and maximum dust-to- 
gas ratios for St = 1 (top) and St = 0.05 (bottom) particles. The color 
represent the different initial dust-to-gas ratios: 1:100 (black), 1 :1000 
(red), and 1 : 10000 (green). More individual super-particles are cap¬ 
tured in the vortices for a low initial e, because their back-reaction is 
less efficient. Since each of these super-particles is less massive than 
with higher £, the overall dust-to-gas ratio for low initial £ is lower 
than that of larger initial £. 

Fig. 6 of the electronic version of this paper). 

The destruction of the vortex occurs on very short 
time scales set by the streaming instability, whose 
growth time is a function of the dust to gas ratio and 
Stokes number of the dust that can be as fast as one or¬ 
bital period (Johansen & Youdin 2007). The new vortex 
amplification on the other hand takes hundreds of or¬ 
bits like in the initial growth as it depends on the radial 
stratification of the disk and the thermal relaxation time 
of gas. 

This vortex disruption and reforming process dis¬ 
cussed here is not to be confused with the vortex in¬ 
stability discussed in Chang & Oishi (2010). They an¬ 
alyze the stability of a three dimensional vortex based 
on the density contrast between the interior of a vor¬ 
tex and the ambient medium and determine that if this 
contrast is higher than a few 10% then the vortex will be¬ 
come unstable. Although we see such high dust density 
concentrations inside the vortex they are very localized 
and therefore our vortices remain stable. For them to be¬ 
come entirely unstable according to Chang and Oishi's 
analysis the density contrast needs to be uniformly high, 
spread out over the entire vortex, which is never the 
case in our simulations. 

Figure 7 shows the maximum dust-to-gas ratio of our 
different simulations with different St. We clearly see 
that St = 1 particles (dash-dotted green line) have the 
highest concentration. As the particle size decreases, the 
concentration also decreases. It is important to note that 
St = 20 particles do not accumulate inside the vortex. 
The highest dust-to gas ratio is reached for St = 1 parti- 


















Fig. 11.— Vorticity co z in units of the local Keplerian frequency Qq of the gas flow (left), dust-to-gas-ratio (middle) and particle collisional velocity 
(right) for St = 1 particles, £o = 10 -2 and after 200 local orbits. Collisional velocities are higher where high particle concentrations are located, but 
hardly exceed 10 _3 c s . 


cles. Larger and lower St reach lower dust-to-gas ratios. 

We now turn our attention to the spatial distribution 
of the dust concentration. By clustering we understand 
what fraction of the dust takes part in the high over¬ 
densities which are the particles triggering the stream¬ 
ing instability. In Fig. 8 we show what fraction of the 
entire dust content has a specific dust-to-gas ratio. 
The dashed vertical lines indicate the initial dust-to-gas 
ratio so = 0-01 and s = 1. For e.g. St = 1 particles 83.74% 
have accumulated in regions with e > 1 whereas for 
St = 0.05 particles only 2.13% of the entire dust mass 
is concentrated in areas with e > 1. The remaining dust 
is spread out thinner. All results can be see in the 7th 
and 8th column of Table 1. Note that for the cluster¬ 
ing we considered a temporal average, but for the max¬ 
imum dust-to-gas ratio £ max (5th column in Table 1) we 
took the absolute maximum form the entire run. There¬ 
fore for runs FI and DG1 £ max > 1 at a given moment 
while on a temporal average this value is not reached. 
A dust-to-gas ratio of e = 1 is significant, because that 
is when dust back-reaction on the gas becomes impor¬ 
tant. This means that for St = 1 back-reaction is a re¬ 


quirement if we want to model particle behavior realisti¬ 
cally. However, for larger St it seems that back-reaction, 
and with that the streaming instability, contributes lit¬ 
tle to the overall dynamical behavior of particles. Yet, 
we already saw that there is a significant difference be¬ 
tween the simulations with and without particle feed¬ 
back for St = 0.05 particles (see Figs. 1 and 7). Without 
feedback all particles were accumulated in the vortex 
whereas the maximum dust-to-gas ratio with feedback 
was around 2-8, two orders of magnitude lower.This 
is confirmed by the Fig. 9, where we compare the clus¬ 
tering of St=0.05 particles with (dashed line) and with¬ 
out feedback (solid line). The particle accumulation is 
similar for low dust-to-gas ratios. However, as e in¬ 
creases, streaming sets in for the simulation with feed¬ 
back and thus regulates the overdensities. Without feed¬ 
back the overdensities can grow unhinderedly which re¬ 
sults in about 80% of all particles accumulating in e > 1. 
We want to stress that although the dust-to-gas ratios 
are different depending on whether back-reaction is in¬ 
cluded or not, it might be difficult to distinguish such 
extremely large local confined concentrations observa- 


















Fig. 12.— Vorticity co z in units of the local Keplerian frequency Qq of the gas flow (left), dust-to-gas-ratio (middle) and particle collisional velocity 
(right) for St = 0.05 particles, £o — 10 -4 and after 200 local orbits. Collisional velocities are higher where high particle concentrations are located, 
but hardly exceed 10 -4 c s . 


tionally with e.g. ALMA (Lyra & Lin 2013). Also for 
St = 0.01 particles that only reach £ max ~ 1 we already 
saw that there is an effect on the gas (see Fig. 5). In this 
context one has to stress that growth rates and unstable 
wavelengths towards the streaming instability depend 
on the St as well as on the actual dust-to-gas ratio. A de¬ 
tailed study of the non linear evolution of the streaming 
instability and the clumping caused by the streaming in¬ 
stability can be found in Johansen et al. (2009) who can 
afford a resolution 10 times higher than we do. One can 
interpret the simulations of Johansen et al. (2009) as a 
detailed study on the dust behavior inside a vortex con¬ 
centration and therefore our vortices generate the dust- 
to-gas overdensities needed for the Johansen et al. (2009) 
scenario. 

We conclude that as soon as even a fraction of particles 
approaches e « 1 back-reaction from the particles onto 
the gas needs to be included to accurately model their 
behavior. 


4.3. Lower Dust-to-Gas ratios 


So far the considered initial dust-to-gas ratio always 
was £o = 1 • 100/ but it has been shown that plan¬ 
ets can also form in low metallicity disks (Mordasini 
et al. 2012; Niedzielski et al. 2009). It is also possible 
that planetesimal formation is still ongoing in a disk 
where some planets have already formed like in the 
51 Peg system (Dumusque et al. 2012) and therefore 
the disk can be depleted in solids. A third reason to 
consider lower initial dust-to-gas ratios is the possibil¬ 
ity that only a small fraction of the entire dust con¬ 
tent is in the particle size regime that we study as a re¬ 
sult of the coagulation-fragmentation balance as stud¬ 
ied by Birnstiel et al. (2012). To account for this we 
perform simulations with lower initial dust-to-gas ra¬ 
tios like £o = 1 : 1000 and £o = 1 : 10000. The number of 
super-particles we put into the domain stays the same, 
while each super-particle represents less mass than in 
cases with higher initial dust-to-gas ratio. 

Figure 10 shows the particle accumulation and result¬ 
ing dust-to-gas ratios for St = 0.05 (simulations F2, DG1, 
DG2) and St = 1 (simulations F3, DG3, DG4). In cases 
with low £q more super-particles concentrate in one lo- 
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Fig. 13.— Gas vorticity co z in units of the local Keplerian frequency 
i?o and dust-to-gas ratio for St=l particles with twice the standard 
resolution. Streaming is still clearly visible and there is no significant 
difference with respect to the lower resolution simulation (see Fig. 3). 


cation. Since each of these super-particles has less mass 
compared to in simulations with higher dust-to-gas ra¬ 
tio, the back-reaction is less effective. Thus the absolute 
dust-to-gas ratio values to be reached in low metallicity 
systems are almost as high as in high metallicity sys¬ 
tems. The gas is not affected as much by the dust par¬ 
ticles as in previous simulations. Thus the vortices, al¬ 
though still disrupted slightly by the back-reaction, are 
no longer torn apart. The particles are trapped more 
tightly and cannot leave the vortex. 

We conclude that the relative concentration in low 
metallicity disks is much stronger than in high metal¬ 
licity systems, thus the absolute dust-to-gas ratio val¬ 
ues to be reached in low metallicity systems is almost as 
high as in high metallicity systems. In the end always 
a roughly identical dust-to-gas ratio is reached (within 
one order of magnitude edg), e.g. at the physical condi¬ 
tion for triggering the streaming instability. 

5. COLLISIONAL VELOCITIES 

In our simulations collisions between the particles are 
neglected. Already for the physical conditions of a pro¬ 
toplanetary disk with real particle numbers a typical 



Fig. 14.— Dust-to-gas ratio for St=l particles, an initial dust-to-gas 
ratio of £o = 3 and no vortex. Clumping of particles still happens, 
therefore the streaming instability is not an effect of the vortex captur¬ 
ing alone. It can also be resolved with standard resolution and without 
radial stratification. 


collision time, i.e. the time after all particles have col¬ 
lided once is typically a few thousand orbits. With our 
super-particle approach we have no chance to ever ob¬ 
serve a real collision in the course of our simulation. Still 
collision speeds are important to know to study possible 
growth or disruption in statistical codes (e.g. Birnstiel 
et al. 2012). 

So when we talk about collisional velocities we rather 
mean relative velocities of neighboring particles from 
which we extrapolate likely collision speeds. In the 
following we will discuss how we calculate their colli¬ 
sional velocities. 

Since we use a two-dimensional approach, our model 
can only provide first estimates for collisional velocities. 
For instance we neglect small scale turbulence that will 
arise from the vertical structure of the disk and particle 
layer (e.g. Kelvin-Helmholz: see Johansen et al. 2006a) 
and elliptic instability due to 3D waves propagating in¬ 
side the vortex, but also from the tail of the Kolmogor- 
rov cascade that is not present in 2D simulations. In 
general 2D flows tend to create large vortical structures 
rather than allowing energy to cascade down to the dis- 































14 


Raettig et al. 


sipation scale. 

Therefore we take a rather simple approach to esti¬ 
mate collisional velocities, which are basically the local 
dispersion in the particle velocities. More exactly: We 
consider each grid cell separately, since only particles 
which are close together have a chance to collide. From 
each individual particle's velocity in a given grid cell 
the mean particle velocity of all particles within this grid 
cell v is subtracted. Afterwards the average of these 
residual velocities within one grid cell is taken. This 
results in the mean collisional velocity within this grid 
cell. To get the average collisional velocity of the entire 
system v co \\ we average over all grid cells that include 
at least 10 particles. The results can be seen in Table 2 
for St=l and St=0.05 particles. Figs. 11 and 12 show the 
spacial distribution of the collisional velocities. Higher 
collisional velocities coincide with regions of high parti¬ 
cle density. 

For tightly coupled particles (low Stokes numbers) it 
is neccessary to first subtract the local gas velocity at the 
position of a particle from Otherwise we would mea¬ 
sure the velocity dispersion of the gas rather than of the 
particles. This results in low collisional velocities Ai? co n, 
for the smaller St=0.05 particles. This is not so much the 
case for St=l particles as they only marginally couple to 
the gas. 

If we compare these collisional velocities to the the¬ 
oretically predicted collisional velocities by Ormel & 
Cuzzi (2007) for a given turbulent oc value we see grave 
deviations. Already our gas velocity u rms , which is 3 
to 4 times as large as the turbulent gas velocity ut = 
\foLC s , where a: is a dimensionless viscosity parameter 
(Shakura & Sunyaev 1973), suggest a deviation from 
their model. For St = 1 particles they predict a ratio of 
particle collisional velocities to turbulent velocities of 1, 
yet our value is of the order of 10 -3 . This shows that 
the theory of Kolmogorov turbulence is not applicable 
to the vortex environment, at least as long as we neglect 
the third dimension and thus cannot predict particle col¬ 
lisions inside vortices accurately. 

From our model we see that the oc values measured 
outside the vortex are not a good proxy to estimate the 
collisional velocities inside the vortex. In our 2D simu¬ 
lations the vortex may appear calmer than it would ac¬ 
tually be in 3D where elliptical instability can produce 
turbulence inside the vortex core (Lesur & Papaloizou 
2010; Lyra & Klahr 2011). 

Yet, as we stated before, this is a rather simple ap¬ 
proach and these collisional velocities can only be seen 
as a first estimate. There are several factors that can 
increase the velocities. Higher resolution and a treat¬ 
ment of three dimensions will lead to properly resolving 
the streaming instability and with this particle velocities 
will increase (Johansen et al. 2009). Including the verti¬ 
cal stratification will increase particle velocities, too, as 
also the vortex dynamics become more complex. 

It is not possible to quantify the error introduced by 
our two-dimensional approximation. For that purpose 
we have to wait until we have evaluated full three- 
dimensional models, which are numerical wise signif¬ 
icantly more expensive. 


6 . PARTICLES COLLECTED INSIDE THE VORTEX 
AND POSSIBLE GRAVITATIONAL COLLAPSE 

Our 2D simulations tell us that the dust-to-gas ra¬ 
tio is increased significantly inside a vortex. A de¬ 
tailed study on gravitational fragmentation of the par¬ 
ticle layer would require 3D modeling as the criterion 
for collapse is to reach a critical 3D density of Pr, which 
is the Roche density. The Roche density is defined as the 
density of a clump that cannot be sheared apart by tidal 
forces from the central object. 

It can be derived by equating the the gravitational ac¬ 
celeration on the surface of a clump 


where M c and R c are the clumps mass and radius re¬ 
spectively, with the tidal acceleration at on the clumps 
surface. The tidal acceleration is the difference between 
the gravitational pull of the central star at the clumps 
surface and its center of mass. To first order we get 

a t « =f2R c ^L, (15) 

where M* is the mass of the central object and r is the 
distance between star and particle clump. If we now 
equate these two accelerations, exchange the masses by 
their respective densities and solve for the clump den¬ 
sity we get 

Pcrit > PR — 2 nr 3 • 

If the density of the particle clump is higher than this 
critical density it will be held together by its own grav¬ 
ity. If this is not the case then the particle clump will be 
torn apart by tidal forces. 

There are additional forces that act to prevent col¬ 
lapse besides gravitational tides, e.g. erosion and tides 
exerted by the vortex (Lyra et al. 2009b), but we only 
aim for a first estimate here and therefore neglect these 
forces. 

If we take the mass of the Sun for the mass of the 
central object then we get pr(1 AU) = 2.83 x 10 _7 gcm“ 3 
and pr(5.2AU) = 2.01 x 10 -9 gcm -3 . The highest over¬ 
densities we observed reached for St = 1 particles is 
Zd ~ 80Zg. This means that the trapping alone will in 
most cases not lead to Roche densities, yet we have ne¬ 
glected sedimentation. As was shown in Johansen et al. 
(2009) an increase in the metallicity of the disk by a fac¬ 
tor of a few is sufficient to trigger both streaming insta¬ 
bility and in succession gravitational collapse into plan- 
etesimals. 

For St=0.01 particles for example it is apparent from 
Fig. 5 that the material inside the vortex structure 
reaches dust-to-gas ratios of e ~ 0.1. This is a concentra¬ 
tion in local available dust by a factor of 10. As shown 
by Johansen et al. (2009) a concentration by a factor of 3 
would be sufficient to trigger planetesimal formation. 

Nowhere outside the vortex these e values are 
reached. Comparing this to the clustering in Fig. 8 we 
see that about 40% of the entire particle fraction are in 
concentrations of e > 0.1 and thus about 40% of the en¬ 
tire particle fraction is captured inside the vortex. Do¬ 
ing the same analysis for the other particle sizes even 
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leads to values up to 98% which means that the captur¬ 
ing mechanism is very efficient. Note that these are best 
case values, since in our simulation no material is lost. 
If a particle is not captured by the vortex the first time it 
passes it, it can be captured the second time due to the 
periodic boundaries, whereas in a real disk this particle 
would be lost to the vortex. On the other hand we also 
neglect the influx of particles from larger radii (see Birn- 
stiel et al. 2012), thus one can argue that our fixed global 
metallicity is a conservative under estimation, as even 
more material might end up in the vortical trapping re¬ 
gion. 

7. SUMMARY AND CONCLUSION 

In this paper we analyze how particles and vortices 
sustained by the radial stratification of a disk affect each 
other. In particular, we investigated whether particles 
of various sizes can concentrate in vortices in the pres¬ 
ence of particle feedback, and whether vortices remain 
a long-lived phenomenon in this case. 

We have conducted two sets of simulations for a series 
of Stokes numbers (St = 20,1,0.05 and 0.01): one where 
there is only gas drag on the particles and one where 
particles also exert drag on gas. This becomes important 
if the initial dust-to-gas ratio is locally enhanced from 
0.01 to about 1. 

We see that, without back-reaction, particles of St = 
0.05, and St = 1 are swept up very efficiently by the 
vortices. For St=l the concentration is so efficient that 
all particles end up in a single grid cell. Smaller and 
larger particles can escape from the vortex again. For 
smaller particles this is because they are too slow to fol¬ 
low the dynamical evolution of the vortices. Larger par¬ 
ticles leave the vortex in the azimuthal direction, yet get 
trapped in the zonal flow correlated with the radial po¬ 
sition of the vortex. 

The dynamics changes considerably when particle 
feedback is included. St = 20 particles are hardly af¬ 
fected by the vortex structure and no critical parti¬ 
cle concentration to trigger the streaming instability is 
reached. Therefore the gas (and thus the vortex) is 
not affected by these large particles. In conclusion, if 
smaller particles coagulate to sizes corresponding to 
St=20, these bodies would neither be captured by vor¬ 
tices nor would they be retained inside them. Yet we 
highlight that our simulations show these bodies of 
large Stokes number getting trapped in the zonal flow 
related to the vortex. This effect has already been ob¬ 
served in global simulations (Figure 11 in Lyra et al. 
2009b). 

For smaller sizes we see high concentrations inside 
the vortex. Particles of St = 1 concentrate very efficiently 
in the center of the vortex, leading to streaming insta¬ 
bility. The instability leads to strong vorticity perturba¬ 
tions, that eventually disrupt the vortex. With decreas¬ 
ing particle density the vortex can be re-established and 
the cycle repeats. The dust-to-gas ratio can be locally in¬ 
creased up to e « 80. We measure that over 80% of the 
dust is concentrated in grid cells with e > 1. 

Smaller particles are not that strongly concentrated. 
Instead they are spread out over the entire vortex and 
take part in the streaming instability (Youdin & Good¬ 


man 2005; Youdin & Johansen 2007). This instability 
typically shows up when e > 1 even when only a small 
fraction of the available dust takes part in the instabil¬ 
ity. For St = 0.05 particles only 2.13% of the entire mass 
reached e > 1. The maximum dust-to-gas ratio reached 
are e ~ 3 for St = 0.05 and e ~ 1 for St = 0.01. 

Although the particle concentrations achieved are 
quite different for the different particle sizes, the over¬ 
all mass of particles accumulated in a vortex is roughly 
the same which matches quite well the analytical predic¬ 
tion of Lyra & Lin (2013). Around 90% of the entire dust 
content is swept up by the vortex. This corresponds to 
a dust density about four times higher inside the vor¬ 
tex than in the background, and about 30 times higher 
than in the region outside the vortex. This highlights the 
prospects of detecting vortices by their enhanced dust to 
gas ratio (Wolf & Klahr 2002; van der Marel et al. 2013). 

We also conducted simulations with lower initial 
dust-to-gas ratios, for three reasons. These are: 1) to 
see whether planetesimal formation can be triggered 
around low metallicity stars; 2) to see if evolved dust 
populations in disks (with significant dust mass in at 
least planetesimal-mass objects or already accreted due 
to radial drift) can still form planetesimals, and; 3) to 
mimic a situation in which only a small fraction of the 
total dust content grows to a size of St = 0.01 - St = 1 
as a result of the coagulation-bouncing-fragmentation- 
drift balance (Birnstiel et al. 2012). We see that although 
the initial dust-to-gas ratio £o is a factor 10 or 100 lower, 
the locally reached maximal dust-to-gas ratio is still of 
the same order of magnitude. The relative particle con¬ 
centration increases until the back-reaction becomes sig¬ 
nificant, which is independent of the initial dust-to-gas 
ratio. 

The first estimate of the collisional velocities of the 
dust particles are much lower than what is to be ex¬ 
pected based on Kolmogorov turbulence of equivalent 
a. These low velocities are possibly a result of our two- 
dimensional approximation and will have to be tested 
in three-dimensional simulations to accurately calculate 
collision speeds and its effect of limiting particle accu¬ 
mulation. 

We confirm that baroclinic vortices are an impor¬ 
tant mechanism for accumulating particles even when 
feedback via drag is considered. The concentrations 
achieved are, depending on particle size, a factor 100 
to 10000 higher than the average value. Even if there 
is only a low amount of dust present, these high over¬ 
densities can be reached. Streaming instability addi¬ 
tionally enhances the dust concentration. The present 
2D study neglects the vertical structure of the gas disk 
and the dust layer. The expected difference in 3D vs. 
2D simulations lies in the additional concentration ef¬ 
fect by vertical sedimentation that is obviously not pos¬ 
sible in 2D vertically integrated models. Thus, 2D stud¬ 
ies can be said to conservatively underestimate the max¬ 
imum particle concentration that can be reached in 3D. 
On the other hand, in three dimensions, the vortex itself 
is weaker due to elliptic instability (Lesur & Papaloizou 
2009; Lyra & Klahr 2011). Also, due to this instability, 
enhanced particle diffusion is acting inside the vortex. 
These issues will be addressed in a future 3D study. 
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APPENDIX 

TEST OF STREAMING INSTABILITY 

As we stated in Section 4.2 we expect the streaming instability to occur due to the achieved dust densities of 
e > 1. However, with the used resolution we do not resolve all possible modes. We performed two tests in order to 
ascertain that the streaming instability is indeed present. First we simply doubled the numerical resolution, so that 
smaller streaming modes are resolved. For the second test, we modeled a box without a vortex or radial stratification 
(jS = 0), but increased the initial dust-to-gas ratio to a value that would trigger the streaming instability (e = 3). 

There was no significant difference between our standard resolution and the first test case (Figs. 13 and 3, lower 
right). We see the same streaming structures in the particle density, and the achieved overdensities in both cases are 
comparable. Therefore we can deduce that the standard resolution we used for all simulations is high enough to 
capture the essence of the streaming instability. 

The second case also showed streaming structures (Fig. 14). This shows that our resolution is high enough to 
enable the streaming instability, and that the streaming structures are not simply a result of particle accumulation in 
the vortex. The vortex accumulation first produces a high concentration of particles, that is in turn amplified by the 
streaming instability, further increasing the local dust-to-gas ratio. 
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